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D.M. Anderson, G.B. McFadden and A. A. Wheeler* 
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Abstract 

We review the development of diffuse-interface models of hydrodynamics and their 
application to a wide variety of interfacial phenomena. These models have been ap- 
plied successfully to situations in which the physical phenomena of interest have as- 
sociated with them a length scale commensurate with the thickness of the interfacial 
region, (e.g. near-critical interfacial phenomena or small scale flows such as those occur- 
ring near contact lines), and fluid flows involving large interface deformations and/or 
topological changes (e.g. breakup and coalescence events associated with fluid jets, 
droplets, and large-deformation waves). We discuss the issues involved in formulating 
diffuse- interface models for single- component and binary fluids. Recent applications 
and computations using these models are discussed in each case. Further, we address 
issues including sharp-interface analyses that relate these models to the classical free- 
boundary problem, computational approaches to describe interfacial phenomena, and 
models of fully- miscible fluids. 


* Permanent address: Faculty of Mathematical Studies, University of Southampton, Highfleld, Southamp- 
ton S017 1BJ, UK. 
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1. INTRODUCTION 

The nature of the interface between two fluids has been the subject of extensive investigation 
for over two centuries. Young, Laplace and Gauss, in the early part of the 1800’s, considered 
the interface between two fluids to be represented as a surface of zero thickness endowed 
with physical properties such as surface tension. In these investigations, which were based 
on static or mechanical equilibrium arguments, it was assumed that physical quantities, 
such as density, were, in general, discontinuous across the interface. Physical processes such 
as capillarity occurring at the interface were represented by boundary conditions imposed 
there (e.g. Young’s equation for the equilibrium contact angle or the Young-Lap lace equation 
relating the jump in pressure across an interface to surface tension times curvature). Poisson 
(1831), Maxwell (1876) and Gibbs (1876) recognized that the interface actually represented 
a rapid but smooth transition of physical quantities between the bulk fluid values. Gibbs 
introduced the notion of a dividing surface (a ‘surface of discontinuity’) and surface excess 
quantities in order to develop the equilibrium thermodynamics of interfaces. The idea that 
the interface has a non-zero thickness (i.e. it is diffuse) was further developed in detail by 
Rayleigh (1892), and van der Waals (1893) who proposed gradient theories for the interface 
based on thermodynamic principles. In particular, van der Waals gave a theory of the 
interface based on his equation of state and used it to predict the thickness of the interface 
which he showed became infinite as the critical temperature is approached. Korteweg (1901) 
built on these ideas to propose a constitutive law for the capillary stress tensor in terms of the 
density and its spatial gradients. These original ideas have been developed further and refined 
over the past century and we turn the reader to the work of Rowlinson (1979) and Rowlinson 
and Widom (1989), who provide thorough discussions of the historical perspectives and 
complete references to the early work on interfacial and capillary phenomena. 

The notion of a diffuse interface and the use of a capillary stress tensor to model the 
interface between two fluids and the forces associated with it are of central importance to 
the topics under consideration in the present paper. Our focus here is on the use of diffuse- 
interface models which fully couple these notions into a hydrodynamic description. Such 
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models have been used to understand physical and hydrodynamic phenomena that occur 
near a fluid’s critical point. Additionally, developments in modem computing technology 
have stimulated a recent resurgence in the use of the diffuse-interface models for the com- 
putation of flows associated with complex interface morphologies and topological changes, 
such as droplet breakup and coalescence and the highly nonlinear development of classical 
hydrodynamic instabilities. 

In the classical fluid mechanical approach, the interface between two immiscible fluids is 
modeled as a free boundary which evolves in time. The equations of motion which hold in 
each fluid are supplemented by boundary conditions at the free surface which involve the 
physical properties of the interface. This formulation results in a free-boundary problem 
(Lamb 1932, Batchelor 1967, Lighthill 1978, Drazin & Reid 1981, Davis 1983). 

Specifically, in the free-boundary formulation it is assumed that the interface has associ- 
ated with it a surface tension, which on applying a stress balance at the interface gives rise 
to the interfacial boundary condition 

tr-n\±= n/Kn, ( 1 ) 

which relates the jump in the stress across the interface to the interfacial curvature. Here 
a is the stress tensor, n is the unit vector normal to the interface, 7 is the surface tension 
(here assumed to be constant) and 1C is the appropriately signed curvature. In addition, an 
interface between two immiscible fluids is impermeable in which case conservation of mass 
across the interface gives that 

«n|_ = t>n| + = V n , ( 2 ) 

where v represents the velocity of the fluid and V n is the normal velocity of the interface. 
Finally, for viscous fluids, there is continuity of tangential velocity across the interface 

[v — (v • n)n]|* = 0. (3) 

The free-boundary description has been a very successful model in a wide range of sit- 
uations, however there are also important instances where it breaks down. In short, as a 
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physical model it breaks down when the interfacial thickness is comparable to the length scale 
of the phenomena being examined. For example, (i) in a near-critical fluid the thickness of 
the interface diverges at the critical point (Stanley 1971) and consequently the representation 
of the interface as a boundary of zero thickness may no longer be appropriate, (ii) the motion 
of a contact line along a solid surface involves a detailed consideration of the fluid motion in 
the vicinity of the contact line, and may require the treatment of length scales comparable to 
that of the interface thickness, and (iii) the free-boundary description may not be adequate 
for situations involving changes in the topology of the interface (e.g. the breakup of a liquid 
droplet), since these processes fundamentally involve physical mechanisms acting on length 
scales comparable to the interface thickness. In addition to the above situations, another 
difficulty associated with the free-boundary formulation arises in its use in computational 
settings when the free boundary shape becomes complicated or self-intersecting. 

Diffuse-interface models provide an alternative description in the face of these difficulties. 
Quantities which in the free-boundary formulation are localized to the interfacial surface are 
recognized to be distributed throughout the interfacial region. For example, surface tension 
in the classical model is a representation of a distributed stress within the interfacial region. 
In this spirit a continuum theory of the interface may be developed where the reversible part 
of the stress tensor, C, that is associated with surface tension is expressed in its simplest 
form as 

C oc (pV’p + 1 1' Vp| 3 ) I ~ Vp ® Vp, (4) 

where p is the fluid density and the components of the outer product Vp® Vp are given by 
{dpi dxi){dp / dxj ) . The existence of such a stress tensor, often called the capillary tensor, 
was first described by Korteweg (1901). The derivatives of the density that appear in the 
stress tensor arise from the non-local interaction of the molecules within the interface. In 
this situation, the density p is the variable which distinguishes the bulk fluids and the 
intervening interface. In this role, it is known as the order parameter. In contrast, for a 
binary fluid undergoing spinodal decomposition, the composition c naturally plays the role 
of an order parameter (Cahn h Hilliard 1958). Alternatively it is possible that neither the 
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density or the composition may be an appropriate or convenient order parameter; such is 
the case in solidification models (e.g. Caginalp 1985, 1986, Langer 1986). Here it is possible 
to introduce an alternative order parameter, the so-called phase field <f> } to characterize the 
phases. The phase-field assumes distinct constant values in each bulk phase and undergoes 
rapid but smooth variation in the interfacial region. The phase field can be regarded as 
a mathematical device that allows a reformulation of the free-boundary problem and has 
been used successfully in many instances. In particular, phase-field models of solidification 
have been used to compute complicated realistic interfacia 1 structures such as those pret>.;i_t 
during dendritic growth (Kobayashi 1993, Wheeler et al 1993, Warren & Boettinger 1995) 
and Ostwald ripening (Warren & Murray 1996). 

Inherent in diffuse-interface models is an interfacial width which is characterized by the 
length scale over which the order parameter changes. By considering the asymptotic limit 
in which the interfacial width is small compared to the macroscopic length scale associated 
with the motion of the two bulk fluids (i.e. the sharp-interface limit) the diffuse-interface 
model can be related to the free-boundary problem. 

In Section 2 we formulate the diffuse-interface theory of a single-component fluid near its 
critical point. Here we discuss developments for equilibrium and nonequilibrium situations 
and also review the applications addressed with this model. In Section 3 we review the 
developments of these ideas for a binary fluid. In Section 4 we discuss related topics including 
the sharp-interface limit analyses, computational methods for fluid interfaces and models of 
miscible fluids. 

2. A SINGLE-COMPONENT FLUID 

Diffuse-interface models of a single- component fluid have been developed largely from the 
perspective of critical phenomena. While they have been used to study phenomena associated 
with the critical point, they have also been applied to situations away from the critical 
point. With this in mind, we shall review in this section the ideas and applications of 
diffuse-interface models of a single-component fluid. 




DIFFUSE-INTERFACE METHODS IN FLUID MECHANICS April 1995 


Equilibrium 

We begin by first considering the equilibrium state of a non-uniform single- component fluid. 
We assume that an isothermal fluid near its critical point has associated with it the Helmholtz 
free energy functional given by 

7 = j v [pKp,T) + ^K\V P \^dV, (5) 

where V is a control volume, /(p, T) is the bulk free energy density (per unit mass), K is a 
gradient energy coefficient (assumed for simplicity to be constant) and T is the temperature. 
In a simple model the term pf p,T) is assumed to take the form of a double-well with 
respect to the density below the critical temperature and a single-well above the critical 
temperature. The square-gradient term is associated with variations of the density and 
contributes to the free energy excess of the interfacial region, which defines the surface 
energy (Cahn & Hilliard 1958). The form of the energy density can be interpreted in the 
context of statistical mechanics in which the square gradient term arises from attractive 
long-ranged interactions between the molecules of the fluid and in which the gradient energy 
coefficient, if, can be related to the pair correlation function, (e.g. see Irving & Kirkwood 
1950, Bearman & Kirkwood 1958, Yang et al 1976, Bongiomo et al 1976, Abraham 1979, de 
Sobrino & Peternelj 1982 and Davis & Scriven 1982 for further details). 

The equilibrium conditions axe obtained by minimizing T subject to a constraint of 
constant mass, Af, where 

M = f v pdV. (6) 

This leads to the Euler-L* grange equation 

KV*p-(f>f), + \ = 0, (7) 

where A is the Lagrange multiplier associated with conservation of mass. We observe that 
the integrand of T as well as that of the mass constraint are independent of the spatial 
coordinates. Consequently, it follows from Noether's Theorem (Goldstein 1980) that there 
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is a corresponding conservation law given by 


V • T = 0, 


( 8 ) 


where T is a second-rank tensor given by 

T = £/ - v ' , ®aW (9) 

and L = p/(p, T) + |V/>| 2 — A p. Using the Euler-Lagrange equation (7) to eliminate the 

Lagrange multiplier we find that 


-h 


P + KpV i p+-K\Vp\ a 


I — KVp ® Vp, 


( 10 ) 


where p = p 2 f p has been identified as the thermodynamic pressure, see, e.g., Callen (1985). 
Using the divergence theorem, equation (8) may be expressed as 


J T ■ ndA = 0, (11) 

where S is the boundary of a control volume with unit normal vector to. This is the gen- 
eralization to three dimensions of the first integral of the Euler-Lagrange equation (7) in 
one dimension. Equations (8) and (11) suggest that T represents a stress tensor (up to an 
additive divergence-free contribution). We show in the next section that T represents the 
reversible part of the stress tensor. Similar equilibrium conditions have been obtained by 
Blinowski (1973a, 1973b) from the point of view of elastic fluids. A review article covering 
a variety of aspects of this and related theories can be found in Davis & Scriven (1982). It 
was noted by Dunn & Serrin (1985) that consistency with nonequilibrium thermodynamics 
requires a more specific form for the capillary tensor than that used originally by Korteweg 
(1901) and later in the mechanical equilibrium theories of Aifantis & Serrin (1983a, 1983b). 

The equilibrium density profile p(z) obtained using a van der Waals equation of state 
(Callen 1985, Rowlinson and Widom 1989) represents a smooth transition from one bulk 
density to the other (in the two phase region) over a length scale associated with the gradient 
energy coefficient. Such an interface has associated with it a surface energy given by 

7 -'£(£)'*• < 12) 
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Non-equilibrium 


We now pursue the non-equilibrium situation by outlining a thermodynamic procedure in- 
volving local balances of mass, linear momentum, energy and entropy consistent with the 
inclusion of a square-gradient energy term in the internal energy functional. The total mass, 
M, total momentum, V, total internal energy, £, and total entropy, S associated with a 
material volume ft(t) are 

M 


/ prfv, 

(13a) 


[ pvdv ; 

Jm 

(13b) 

J m &«’ + *'■'> + 

(13c) 

Lr™' 

(13d) 


where v is the fluid velocity, e and a are the internal energy and entropy per unit mass, and 
Ke is the gradient (internal) energy coefficient which we assume to be constant. Here we 
have for simplicity neglected body forces such as gravity. The associated physical balance 
laws can be expressed as 

dM 


= 0, 

Lt i 


m • hdA , 


dt 
dV 

dt Jtn(t) 

= f [v • m • n - q E • n] dA, 

dt Jso{t) 

^ + / qs • hdA = f i^dV > 0, 

dt JsQ(t) Ju(t) 


(14a) 
(14b) 
(14c) 
(14d) 

Jm(t) ‘ ->n(0 

where 6Sl(t) is boundary of m is the stress tensor, q E and qs are the internal energy and 
entropy fluxes, respectively, and a proi is the volumetric entropy production. The quantities 
m, qg and qs, which in general include both classical and nonclassical contributions, shall 
be specified below such that their forms guarantee that a vroi is non-negative, as required 
by the Second Law of Thermodynamics. Equation (14a) simply represents conservation of 
mass. Equation (14b) states that the change in total momentum is related to the forces 
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on the boundary (again note that we have neglected body forces). Equation (14c) states 
that the change in total energy is related to the rate of working done by the forces on the 
boundary and also the energy flux through the boundary. Equation (14d) states that the 
change in total entropy plus the entropy flux through the boundary must be equal to the 
entropy production. The definitions (13) and balance laws (14) may be manipulated to show 
that 

r* - (m ~r >:W + (* + *•&') ■ v (?) + v ■ (* - * - ¥&') • (15) 

where T is the reversible part of the stress tensor given by equation (10) with K replaced by 
Kg, and we have used the thermodynamic relationship de = Tds + ( pi p 2 )dp . The following 
specifications ensure that the entropy production is positive 


m = T + r, 

(16a) 

cC 

> 

Sis 

£ 

i 

> 

I 

II 

(16b) 

„ jbvr 

(16c) 

qs T 1 


where k is the thermal conductivity and r is the viscous stress tensor given in the standard 
m ann er as r = ij( V • v) + /i(Vv+ Vu 2 *), where 7} and p are coefficients of viscosity (e.g. Batch- 
elor 1967). This prescription for m is similar to that postulated by Korteweg (1901). In 
this formulation it is assumed that the viscosity and thermal conductivity are, in general, 
functions of the density. We observe that the energy flux qs involves both the classical 
contribution corresponding to the Fourier Law for heat conduction (Carslaw & Jaeger 1959, 
Kittel & Kroemer 1980) and a nonclassical contribution. This nonclassical contribution was 
referred to as ‘interstitial working’ by Dunn & Serrin (1985), who also noted that there is 
no corresponding nonclassical term in the entropy flux. When a square-gradient term is in- 
cluded in the definition for total entropy (13d), a nonclassical entropy flux arises (e.g. Wang 
et al 1993, Wheeler et al 1996, Anderson & McFadden 1996). Using the above forms (16) 
for the stress tensor and fluxes the local balance laws may then be written as 

= -pVv, (17a) 


Dp 

Dt 
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= Vm, (17b) 

ptg. = V • (kVT) + {—pi + r) : Vt>, (17c) 

= V • (JfcVT) + r : W, (17d) 

These equations describe viscous, compressible, nonisothermal flow. In order to solve these 
it is necessary to supply an equation of state. In the isothermal situation, for example, one 
might specify the pressure p through a van der Waals equation of state. 

This and similar models have been developed and studied by a number of authors from 
a variety of perspectives. Fixman (1967) developed a diffuse- interface hydrodynamic model 
in which the reversible part of the stress tensor was identified using a mechanical principle. 
In the model of Felderhof (1970), a more general form of the stress tensor was obtained 
which was compatible with the underlying Lagrangian. Langer & Turski (1973) derived 
a diffuse-interface model that they showed, through a ‘coarse-graining* argument, could 
be related to a molecular model. The model employed by Jasnow & Vinals (1996) was 
derived using a Hamiltonian description. Jacqmin (1996) described a model which also 
included a wall potential to model the interaction between the fluid and a solid boundary. 
Truskinovsky (1993) derived a similar model which also included an additional non-conserved 
order parameter and its gradients. Antanovskii (1996) presented a derivation of the model 
based on a maximum entropy principle. The derivation outlined above is most similar to 
that described by de Sobrino (1976), Dunn & Serrin (1985), Dunn (1986) and Anderson 
& McFadden (1996, 1997). The work of de Sobrino begins in a more general framework 
and invokes symmetry and invariance principles to simplify the gradient dependence of the 
stress tensor to that which is described above. Dunn & Serrin’s approach is similar, but 
mathematically more rigorous and is given from the viewport of rational mechanics. 

A detailed numerical analysis of a simplified version of the diffuse-interface model has 
been performed by Affouf and Caflisch (1991), who present numerical solutions representing 
phase transitions, shocks, and rarefaction waves connecting far-field states and analyse their 
stability. 
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Applications 

CRITICAL POINT SCALING LAWS Extensive analysis using renormalization-group 
techniques have been performed on a diffuse-interface model, commonly known in the liter- 
ature as ‘Model IP (Hohenberg k Halperin 1977), which describes the dynamics of a binary 
fluid phase transition as well as a single- component fluid near its critical point (e.g. Halperin 
et al 1974, Siggia et al 1976, and Hohenberg k Halperin 1977). Such analyses have identified 
diver: *nt transport coefficients and scaling relations associated with near-critical fluids. 

SHEAR FLOWS IN NEAR-CRITICAL FLUIDS Onuki k Kawasaki (1979) and Onuki 
et al (1981) have studied the dynamics of a near-critical fluid in a shear flow using a model 
developed by Kawasaki (1970). They investigated the regime in which the equilibrium cor- 
relation length (i.e. the interfacial thickness) exceeds the length scale associated with the 
shear flow. Among their findings is that the critical fluctuations of classical fluids can be 
drastically altered (e.g. they can become highly anisotropic) by shear flows. These ideas 
have also been applied to polymers under shear flows (Helfand k Fredrickson 1989, Onuki 
1989). 

CAPILLARY WAVES The diffuse-interface model has been used to study capillary 
waves (Felderhof 1970, Turski k Langer 1980, de Sobrino k Petemelj 1985 and de Sobrino 
1985). These authors begin with a diffuse planar interface in equilibrium. They discuss 
capiiir/y waves by means of a linearized theory about the equilibrium state. The linearized 
governing equations are examined in the long wavelength limit and the dispersion relation for 
capillary waves obtained from the classical ffee-boundary problem (e.g. Landau k Lifshitz 
1959) is recovered. These authors recognized the importance of an isentropic treatment 
to the existence of capillary wsves (in the context of a diffuse-interface model) although 
Felderhof derived the result in che isothermal case as well under somewhat more restrictive 
conditions. 

MOVING CONTACT LINES Seppecher (1996) established the governing equations 
for an isothermal viscous flow near a moving contact line a planar solid wall using a 
diffuse-interface model. In particular, a uniform distribution of ‘double forces 1 was used to 
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describe the interaction between the fluid and the wall. In this analysis, the contact-line 
problem was separated into three regions: an external (outer) region far from the contact 
line where the classical theory of capillarity applies; an inner region very near the contact line 
whose dimensions are so small that the thickness of the interface cannot be neglected; and an 
intermediate region which matches the inner and outer regions. The flow in the inner region 
is compressible while the flows in the intermediate and outer regions are incompressible. 
Examples of the density field and flow field in the inner region are shown in Figs. 1 and 2. A 
key result of this paper was that the force singularity present in classical continuum models of 
moving contact lines (e.g. Huh & Scriven 1971, Dussan V. & Davis 1974, Dussan V. 1979) is 
not present when the interface is modeled as diffuse. This is attributed to mass transfer across 
the interfacial region. Additionally, numerical computations reveal a roughly linear increase 
in the apparent contact angle with the contact-line velocity (no contact-line hysteresis was 
considered), which agrees with the general trends observed experimentally (see Dussan V. 
1979). The moving contact-line problem has also recently been treated computationally by 
Jacqmin (1996) and is described in the next section. 

INTERNAL WAVES IN A NEAR CRITICAL FLUID Anderson h Mcladden (1997) 
have recently employed the diffuse-interface model to describe internal waves in a near- 
critical fluid. These internal waves are present in small (centimeter-scale) containers owing 
to the large compressibility of the fluid near the critical point and have been observed 
experimentally in near-critical Xenon by Berg et al (1996). In the experimental work of Berg 
et al, internal gravity wave frequencies were measured both above the critical temperature 
where a single phase exists and below the critical temperature where two phases, separated 
by an interface, exist. The theoretical development of Berg et al consisted of two separate 
models; one above the critical temperature and another below the critical temperature. In 
contrast to these classical hydrodynamic models, the diffuse-interface approach employed by 
Anderson & McFadden allows for a single model to be applied both above and below the 
critical temperature. In their diffuse-interface model, they use a van der Waals equation 
of state to obtain a base-state density profile. Upon this base-state they introduce linear 
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perturbations to calculate internal wave frequencies for temperatures both above and below 
the critical temperature. Predictions of the internal wave frequency from the diffuse-interface 
model are compared favorably with experimented data and theoretical results of Berg et al. 

DROPLETS AND NUCLEATION Several ? ithors have investigated the nucleation 
of droplets (Blinowski 1974 and DellTsola et al 1995, 1996). The focus of this work was to 
ascertain the effects of interfaci&l thickness on the eation conditions of a droplet since, at 
nucleation, the droplet radius may be comparable to the interfacial thickness. Under these 
circumstances, the classical Laplace-Gibbs theory for the equilibrium radius of a droplet is 
called into question. DellTsola et al used an equilibrium formulation of a diffuse- interface 
model developed in earlier work (DellTsola & Kosinski 1993) to study nucleation of spherical 
droplets. In particular, they noted that for microscopic droplets the difference in mechanical 
pressures inside and outside the droplet (which is what is measured experimentally) is not 
the same as the difference between the thermodynamic pressures inside and outside the 
droplet. The mechanical pressure involves a (stress) contribution from the spatial density 
variation (e.g. the term in (10) involving pV 7 p) which at the center of the microscopic bubble 
is important. Based on these ideas, they carefully define quantities such as surface tension 
and droplet radius in a way which generalize the classical notions. Relationships between 
surface tension, droplet radius and the critical nucleation radius are obtained using a number 
of different equations of state. Results for the minimal nucleation radius are compared with 
experimental measurements. 

INSTABILITIES OF PLANAR JETS Nadiga & Zaleski (1996) studied in the instability 
of a planar jet of viscous, compressible, isothermal liquid issuing into its surrounding gas 
phase. They use a van der Waals equation of state to characterize the system. Their 
calculations focus on the high Reynolds number regime (Re = 800) and they investigate the 
effect of surface tension on the stabilization of the jet. 

SPINODAL DECOMPOSITION IN A PURE FLUID Nadiga & Zaleski (1996) also 
considered the spinodal decomposition of a single-component fluid rapidly quenched from a 
temperature above its critical point to a temperature below the critical point. These compu- 
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tations are isothermal and two-dimensional. They find that the inclusion of hydrodynamic 
effects in the model leads to a predicted growth rate which is slightly enhanced from the case 
when the growth is limited by diffusion. Numerical computations show the domain growth 
and plots of domain size versus time are presented. 


3. A BINARY FLUID 


Non-equilibrium 

We now consider the situation of a binary fluid which consists of two components A and 
B . We denote the composition of component A, expressed as a mass fraction, by c. In this 
setting the composition plays the role of an order parameter which distinguishes the different 
phases of the fluid and in this way is analogous to the density in the single- component fluid 
models. As we shall see shortly, however, there are a number of subtle differences between 
the single-component and binary fluid models. 

The governing equations for the flow of a nonisothermal, compressible binary fluid may 
be developed in a similar manner to those for the single- component case discussed above but 
involve square-gradient contributions in the internal energy functional from the composition 
rather than from the density. This procedure gives that the entropy production is 

= (< + (ft + *«%Vc) • v (1) - ic ■ v (fe) 

+ v(«5-®-^F§jVc +*«,), 
where Tc is the reversible part of the stress tensor given by 

T c = (-P + ^alVcl 2 ) / - K E Vc 0 Vc, 
qc is the mass flux of component A and /x c is the generalized chemical potential given by 

( 20 ) 


(18) 


(19) 


de Kb —a 

Me = q vc. 

dc p 


The specification of m, qg, qs and qc to ensure positive entropy production follows in 
an analogous way to the single-component case. In particular, qc = —D J 7(p e /T) t where 
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D is positive and represents the diffusion coefficient. Consequently, the equation for the 
composition is given hy 

*5 - v 'M? (£■**•)]}■ (a > 

This is well-known Cahn-Hilliard equation, used to model spinodal decomposition by Calm 
(1961) and Hilliard (1970), modified to account for fluid motion. The equations governing 
density, velocity and temperature (energy) are similar to those for a single-component fluid 
given above by equations (17). Two subtle differences between the binary fluid equations 
and the single-component fluid equations are worth mentioning. First, the reversible part of 
the stress tensor T c does not contain the counterpart of the Laplacian term appearing in T, 
see equation (10). This is because the order parameter c is given per unit mass in contrast 
to the density p. Second, the order parameter c is governed by the modified Cahn-Hilliard 
equation (21), whereas the density p is governed by the continuity equation (17a). 

Details of the derivation of the governing equations for a binary fluid have been dis- 
cussed by a number of authors. Blinowski (1975) considered the binary fluid case and also 
more general multi-component systems. Starovoitov (1994) derived a bii^y fluid diffuse- 
interface model using a virtual power method. Antanovskii (1995) derived the equations for 
nonisothermal, viscous, quasi- compressible flow. His derivation was based on a maximum 
entropy principle, similar to the approach described here, but used virtual work arguments 
to identify the forms of the stress tensor and fluxes. Gurtin et al (1996) derived a model for 
an isothermal, incompressible flow using microforce balance laws. Jasnow & Vinals (1996) 
illustrated the derivation of the equations for a binary fluid using a Hamiltonian formalism. 
Although their derivation was given for an isothermal, inviscid, incompressible fluid, the 
model was extended to the case where the fluid was viscous and the temperature field varied 
slowly in time. A similar account can be found in Jacqmin (1996), who discussed the model 
in a ‘potential form’ as well as making note of its relation to the above ‘stress form’. The 
merits of these two equivalent formulations were discussed in terms of specific applications. 
Lowengrub & Truskinovsky (1997) presented a thorough derivation of the diffuse-interface 
model based on entropy production and paid particular attention to the difference between 
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compressible and quasi-incompressible fluids. The quasi-incompressible situation describes 
the case where the fluid density is independent of the pressure. However, the bulk states 
may have different densities and the flow in the interfacial region is in general non-solenoidal 
(V - v ^ 0), resulting in an expansion or contraction flow upon phase transformation. These 
authors argued that within the context of quasi-incompressibility, the appropriate thermody- 
namic description is in terms of a Gibbs free energy, in which the pressure is an independent 
variable determined by the transport equations rather than a quantity determined thermo 
dynamically. 

Applications 

THERMO CAPILLARY FLOWS Antanovskii (1995) used his nonisothermal binary fluid 
model to compute one-dimensional thermocapillary flow in a gap. In this situation, two fluid 
phases, characterized by their different compositions, were separated by a planar diffuse 
interface along which a temperature gradient was imposed. Calculations for different values 
of interfacial layer thickness and viscosity ratios were presented. Jasnow Sc Vinals (1996) 
investigated thermocapillary migration of droplets of one phase in the stirrounding phase. 
They focused on drops with radii on the order of ten times the correlation length. Their 
calculations show the motion of the droplet through the temperature gradient as a function 
of time and the dependence of its velocity on this temperature gradient. Also shown is a 
sequence in which two droplets coalesce. 

SPINODAL DECOMPOSITION Gurtin et al (1996) considered spinodal decomposition 
occurring in an isothermal binary fluid. Their computations, which begin with an initial 
random distribution of the composition, show the coarsening process explicitly. They note 
that the main effect of the hydrodynamic interactions on this process is the flow-induced 
coalescence of droplets. Their calculations specifically show the flow associated with these 
coalescence phenomena. They compute the domain (structure) size, compare with classical 
predictions, and find that for long times the growth is faster than the classical scaling law 
for coarsening by purely diffusive mechanisms. Jasnow Sc Vinals (1996) consider a similar 
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situation but include the effects of a nonunifonn temperature field. Their results, displayed 
in the time-sequence in Fig. 3, show the spinodal decomposition of a binary fluid in a 
rectangular cell across which a vertical temperature gradient is imposed (the cell is hottest 
at the bottom). Here, an initially random composition distribution (top frame) undergoes 
coarsening and domain growth. At an intermediate time (middle frame) and a later time 
(bottom frame), their calculations show nonuniform coarsening wherein the smaller (larger) 
scales occur in the warmer (cooler) regions. 

MIXING AND INTERFACIAL STRETCHING Chella & Vinals (1996) compute the 
mixing and interfacial stretching of an isothermal, incompressible binary fluid in a shear 
flow with equal densities and viscosities in the two phases. An initially planar interface is 
shown to distort under an imposed shear flow. The wrapping up of the two phases is shown 
for different values of the capillary number (surface energy). The amount of interfacial 
stretching is computed and compared with an analytical solution for that associated with a 
passive scalar. An increase in the surface energy corresponds to a decrease in the amount 
of stretching and in this way opposes the effect of the shear flow. Further calculations by 
Chella and Vinals of flow in a driven cavity are shown in Fig. 4. A number of breakup and 
coalescence events can be seen in this sequence. 

DROPLET BREAKUP Jacqmin (1996) calculated the breakup of an inviscid fluid in 
the context of a two-dimensional isothermal model. The initially elongated droplet relaxes, 
oscillates and breaks apart into two separate droplets. In these calculations, the diffusion 
coefficient (D in our notation) was made velocity dependent. In addition to plots showing 
the actual droplet breakup process, the kinetic and surface energy evolutions are also given. 
In the case where the flow is strongly damped by using a larger value of Z), the energy of 
the droplet decreases rapidly and breakup does not occur. 

WAVE-BREAKING AND SLOSHING Jacqmin (1996) also applied the diffuse-interface 
model to a large-deformation sloshing flow in a two-dimensional rectangular domain. Inter- 
action with the solid wall of the container is modeled using a wall potential associated with 
a 90° contact angle between the wall and the interface. The two-phase fluid in a uniform 
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(vertical) gravitational field was subjected to an oscillating horizontal acceleration whose 
amplitude increased linearly in time. The motion of thj initially planar interface between 
the two stably stratified phases is computed once the horizontal accelerations begin. It is 
shown that the shape of interface becomes highly nonlinear and several breakup and coa- 
lescence events can be seen to occur. It is noted that the prediction of coalescence events 
may occur more quickly than is physically realistic since a high level of accuracy is needed 
to resolve the draining layer between the two coalescing phases. 

MOVING CONTACT LINES Jacqmin (1996) investigated the fluid motion near a 
moving contact line. These steady calculations show the existence of a dynamic contact 
angle associated with the interface away from the immediate vicinity of the contact line. 
Abo observed is a streamline (similar to that observed in the experiments by Duss n V. &; 
Davis (1974)) issuing from the contact line region into the displaced fluid. 

NTJCLEATION Lowengrub & Truskinovsky (1997) considered nucleation (and annihi- 
lation) of an isothermal, spherically symmetric equilibrium droplet. In particular, they used 
the diffuse-interface model to describe the situation where the droplet size was comparable 
to the interface thickness. In the incompressible case, where the density was uniform ev- 
erywhere, they performed an analysis of the spherically symmetric Cahn-Hilliard equation 
with a free-energy density composed of piecewise parabolas to find new analytic solutions for 
the compositions. In the compressible case they proceeded numerically and concluded from 
their results that compressibility has little effect on the interfacial structure of the droplet. 

4. RELATED TOPICS 

Sharp-Interface Limit 

As noted in Section 2, the diffuse-interface models may be applied away from the critical 
point, where the interfacial thickness approaches that of a ‘sharp’ boundary. The use of 
the diffuse-interface models in these regimes may be justified by demonstrating that they 
approach asymptotically the free-boundary formulation. This is achieved by adopting what 
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is commonly known as the ‘sharp- interface limit 1 . 

There has been a significant effort in the context of phase- field models of solidification 
to address the sharp-interface limit and to compare the model to well-known results of the 
free-boundary problem (e.g. C agin alp 1989, Braun et al 1994). 

One ov the key features in the diffuse-interface models described here, which is not present 
in the phase-field solidification models studied to date, is the (vector) momentum equation 
which involves the distributed capillary stresses. Consequently, much of the emphasis in 
terms of the sharp-interface analyses of the diffuse-interface model has been to recover the 
classical interfacial boundary conditions associated with the stress balance at the interface. 
Here we describe the efforts in this area and outline a simple reduction of the momentum 
balance to the interfacial stress jump using a pillbox argument. 

Antaaovskii (1995) derived from his diffuse-interface model of a binary fluid the special 
cases cf the classical hydrostatic balance for a flat interface in equilibrium and the Young- 
Laplace equation for a spherical interface in equilibrium. The latter case has also been 
considered for a single-component model by Blinowski (1979). Nadiga and Zaleski (1996) 
also confirmed numerically that their diffuse-interface model accurately recovered the clas- 
sical results for a flat interface and for a liquid droplet in equilibrium. Jasnow and Vinals 
derived from the capillary term in their momentum equation the appropriate sharp-interface 
tangential and normal forces when the surface tension was a slowly varying function along the 
interface. More detailed analytical approaches have been addressed by Starovoitov (1994), 
Anderson & McFadden (1996) and Lowengrub & Ttuskinovsky (1997). 

To illustrate these ideas, we apply a pillbox argument to the momentum balance (17b) to 
show how the classical stress balance at a fluid-fluid interface can be derived from the diffuse- 
interface model. We define a small parameter e measuring the thickness of the interface by 
Ke = We then consider the surface Si defined by the contour of density upon which the 
interfacial region collapses in the limit e — ► 0 and define a pillbox enclosing a portion of this 
surface at a fixed point in time in such a way that the top of the pillbox is above the surface 
at a height r = 6 and the bottom of the pillbox is below the surface at a height r = — 5. 
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Here, r is a local coordinate normal to the interface. The key limit in the pillbox argument 
is that e S <C L where L is an 0(1) length scale associated with the outer flow. In this 
limit, the volume of the pillbox becomes negligible on the outer scales but the variations in 
the density, which define the interfacial region, occur over a region fully contained within the 
pillbox. Also in this limit, the unit normal vectors on the top and bottom of the pillbox are 
ft and —ft respectively while the unit normal on the side is given by rh f note that th « ft = 0). 


We integrate equation (17b) to obtain 


0 = J v ^^ + V {pv®y-V rnjdV t 

(22) 

where we have used (17a) and the fact that 


Dv d(pv) . _ 

p ot = at 

(23) 

Next, we note that 


/ / (pvJui-nsdS, 

JVp Ot 

(24) 


where vj is the velocity of the surface Sj described above and S p denotes the surface of 
the pillbox. This result follows by translating to a frame moving with the interface so that 
d(pv)/dt = d(pv)/dt — vj • V(pv) = d(pv)/dt' — V • (pvj ® v ) + pvV • vj. The terms 
d{pv)jdt + pvV ■ vi are bounded and hence do not contribute as the pillbox volume goes to 
zero. We use this and the divergence theorem to obtain from equation (22) 


0 = J [pv(u — vj) - hs — n$ • mj dS. (25) 

We further argue that the fluid velocity terms are bounded (so that they do not contribute 
to the integral over the side surface of the pillbox) and that the nonclassical terras do not 
contribute to the upper and lower surfaces of the pillbox (since e < f) so that 

0 = J {pv(v — vi) • ft — ft ■ (— pi + r)] + dA - J ^ m • Tcf5, (26) 


where A is the portion of 5/ within the pillbox. Now, local to the interfacial region we have 
Vp ~ p r ft and V a p ~ p TT to leading order. Then m • T ~ (m • T • m)m so that 


0 = j [pt?(v - vj) * ft - ft • (— pi 4* r)] + dA - ^ J (m • T • m)mdrdl } (27) 
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where C is the contour defined by the intersection of Si and the pillbox surface and dl is the 
increment of arclength along C. We next define the scalar 

7 = r {m-T m)dr, (28) 

J—OO 

which can be shown to be equal to the excess Kramer’s potential, that is, the surface energy 
(Anderson & McFadden 1996). We then apply the surface divergence theorem (Weatherburn 
1925) which allows us to write 

£ 7 mdi = £ [V s7 - 7 ( Vs • n) n] dA, (29) 

where Vs 7 is the surface gradient of 7 and V$*n, the surface divergence of n, is the interfacial 
curvature, K. This gives 

0 = jf { \pv{v - vj) • n - n * (- pi + r)] * - [V 57 - 7 £n] j dA. (30) 

Finally, noting that the area of integration is arbitrary yields 

[pv(v - u/) • n - n ■ (-pi + r)]^ = V 57 - 7 ^n, (31) 

which is the classical stress balance at a fluid-fluid interface (Delhaye 1974). Note that the 
first term on the left-hand side represents a jump in the momentum across the interface and 
is zero when there is no mass flux across the interface. 

Other Computational Methods 

The diffuse-interface models share common features with a number of methods developed 
from a more computational point of view. We outline several of these methods below and 
highlight some of their applications. 

The volume of fluid (VOF) method (Hirt & Nichols 1981, Hyman 1984, Tsai & Yue 
1996) is a numerical approach to the free-boundary problem in which an auxiliary function 
F } which distinguishes one fluid from another, is introduced in order to identify the shape 
and evolution of the free boundary. This function satisfies DF/Dt = Oj that is, it is advected 
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with the flow. Each computational cell has associated with it a value of F and those cells 
which take on a value between two bulk values of F are assumed to contain part of the 
interface. The normal to the interface in the cell is determined by the direction of the 
largest local gradient and the position of the interface in the cell is arranged so that F is the 
fractional volume of fluid in the cell. The VOF method, like the diffuse-interface approach, 
allows a straightforward description of flows involving complicated boundary shapes and 
topological changes. In contrast to the diffuse-interface approach, however, free surface 
boundary conditions in the VOP method must still be applied at the free boundary. 

Brackbill et al (1992) developed a ‘continuum surface force* model wherein they iden- 
tified a volume force which represents surface tension spread over a small but finite three- 
dimensional interfacial domain. This volume force was related to a ‘color’ function which, for 
example, can represent density for incompressible flows. The defining characteristics of this 
volume force are that it gives the correct surface force in the limit of a sharp interface and 
is nonzero only in the interfacial region. This approach allows a single-domain description 
of the two-fluid system and does not require the direct application of boundary conditions, 
which are built into the governing equations. A number of numerical results are presented 
for both static and dynamic situations: ( i ) a static drop (rod), («) a nonequilibrium rod 
upon which capillary waves move along the surface, (in) the Rayleigh-Taylor instability and 
the associated interfacial deformation, (iv) flow induced by wall-adhesion whereby the fluid 
conforms to an imposed equilibrium contact angle on the wall, and (v) jet-induced tank 
mixing and liquid reorientation in microgravity environments. 

Another related approach has been developed by Unverdi and Tryggvason (1992a, 1992b). 
Their approach is a front-tracking technique which employs a numerically-diffuse description 
of the interface. They construct an indicator function, based on the known position of the 
(sharp) interface, which identifies fluid properties such as density and viscosity. This function 
is then artificially spread out over a small region on the scale of the computational mesh size, 
allowing the fluid properties to vary smoothly through this interfacial region. The surface 
force (i.e. surface tension) is also distributed over this interfaci&l region so that a single- 
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domain approach can be used to calculate the flow. This flow then determines how the 
interface is advected. In Unverdi and Tryggvason (1992a), both two- and three-dimensional 
multiple bubble motion and interaction are presented. In Unverdi and Tryggvason (1992b), 
the Rayleigh-Taylor instability in two and three dimensions and also bubble-bubble interac- 
tion in three dimensions are computed. Nobari et al (1996) have recently used this approach 
to compute head-on collision of two viscous liquid droplets with surface tension. Here, rup- 
ture is modeled by artificially removing the thin film between the two drops at a prescribed 
time. They found that if no rupture takes place, the drops always rebound, but that when 
rupture occurs the drops may later split. This method has been extended by Juric and Tryg- 
gvason (1995, 1996) to describe flows in the presence of phase change. Here, they applied 
the model to vapor bubble dynamics and film boiling. 

Another highly successful computational scheme that has been applied to interfacial 
motion is the level set method, see Osher & Sethian (1988) and Sethian (1996). With this 
method, the interface is represented as a level set of a smooth auxiliary function which is 
computationally analogous to the order parameter used in diffuse-interface descriptions. An 
advantage of the level set method is that the interface remains sharp in this formulation, 
which eliminates the need for added numerical resolution in the direction normal to the 
interface. Within the context of fluid mechanics and two-fluid flows, surface tension, for 
example, is represented in the momentum equations as a distributed force through the use 
of a smoothed delta-function (Sussman et al 1994, Chang et al 1996). The momentum 
equations are then used to compute the flow over the whole domain and the level set function 
is advected with the flow. After a normalization procedure, the level set function determines 
the new position of the interface. Mulder et al (1992) applied the level set approach (without 
surface tension) to study the Rayeigh-Taylor and Kelvin-Helmholtz instabilities within the 
context of compressible gas dynamics. Sussman et al (1994) used this approach to model 
incompressible two-dimensional rising bubbles and falling drops, which show large distortions 
and topological changes. They also show two impacting droplets as well as a single droplet 
impacting a surface. Chang et al (1996) also applied a level set approach for incompressible 
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fluids to several topologically complex flows. They presented computations of two merging 
fluid bubbles of equal density said also two merging fluid bubbles of different density. Further, 
they investigated the Rayleigh-Taylor instability of an initially motionless, vertical column 
of fluid and show calculations of the subsequent vortex sheet roll-up phenomena. 

Miscible Fluids 

The ideas that are involved in the d ; ffuse-interface models described in the preceding sec- 
tions are similar in many ways to those for miscible fluids. In particular, fluid stresses 
(i.e. Korteweg stresses) that arise due to concentration and density gradients at the interface 
between two miscible fluids lead to the notion of surface tension between miscible fluids. 
Joseph (1990) has investigated these ideas both experimentally and theoretically. He per- 
formed a number of experiments and highlighted other experimental work in which miscible 
liquid droplets rising or falling in another liquid exhibit capillary-type effects; that is, their 
shapes are consistent with the presence of surface tension on the interface. In the theoretical 
development, the equations governing the motion of the fluid are similar to those presented 
in the above section on binary fluids. They are the continuity equation, the Navier-Stokes 
equations modified to account for the gradient stresses, the heat (energy) equation, and a 
standard diffusion equation, rather than the modified Cahn-Hilliard equation (21), that de- 
scribes the evolution of the composition. Besides the difference in the equations governing 
the composition, another key difference between the description given by L^cph (1990) and 
that for a binary fluid described in section 3 is that his model described a generalized (or 
quasi) incompressible fluid. Here it is assumed that the density is a function of compo- 
sition and temperature but is independent of the pressure. Although the density then is 
unchanged by pressure variations, the associated flow may still be nonsolenoidal ( V • v ^ 0) 
due to composition or temperature variations. This notion has been carefully adopted into 
the diffuse-interface model for a binary fluid by Lowengrub & Ttuskinovsky (1997). 

The theoretical model described by Joseph (1990) has been analyzed in more detail by 
Galdi et al (1991). These authors rework the equations and identify a new solenoidal velocity 
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which is a linear combination of the original velocity field and the concentration gradient. 
Associated with thii new velocity is a pressure field that is a linear combination of the 
original pressure and the divergence of the original velocity. They address the linear and 
energy stability of a quiescent, vertically (unstably) stratified incompressible fluid in which 
Korteweg stresses arise due to composition gradients. This is the analog of the classical 
Benard problem with the exception that the authors do not immediately invoke the Boussi- 
nesq approximation. They find that the stability results depend strongly on the value of the 
coefficient in the Korteweg stress (and in particular on its sign). When its sign is chosen 
consistent with that in the capillary stress term described in Section 2, and when its effect 
is strong enough, an unconditionally stable base state is predicted. 

An notable idea described in the paper by Joseph (1990) is that of a dynamic (or nonequi- 
librium) surface tension between two miscible fluids. That is, although for two miscible fluids 
it is not clear that one can define an equilibrium surface tension, based on the idea of stresses 
associated with gradients in density or composition, transient or dynamic surface tension be- 
tween two mixing phases can be studied. Joseph et al (1996) consider this situation in detail. 
In this paper, they carry out an analysis using the model of Joseph (1990) of transient or 
dynamic interfacial tension during the smoothing of an initial discontinuity of composition 
across plane and spherical surfaces separating two miscible liquids. It is found that the 
dynamic interfacial tension decays in time like t~ l S 2 and has contributions from Korteweg 
stresses and from the expansion velocity (i.e. the nonsolenoidal part) which also involves 
the rate of change of viscosity with composition. For a plane mixing front, diffusion has 
a similarity solution and they show that there is no associated pressure jump due to the 
Korteweg terms (i.e. surface tension does not lead to a jump in pressure when the curvature 
of the interface is zero). The only pressure difference across the mixing front is due to the 
hydrostatic pressure difference (i.e. due to gravity) across the layer since its thickness is 
growing with time. For an initially spherical droplet placed in another (miscible) fluid they 
identify a pressure jump associated with two effects: the first is due to the Korteweg stress 
and the second is due to the expansion velocity and is proportional to the rate of change of 


- 24 - 




DIFFUSE-INTERFACE METHODS IN FLUID MECHANICS 


April 1995 


viscosity with composition. The latter term can take on either sign, depending on how the 
viscosity depends on the composition. 

5. SUMMARY 

We have reviewed the development and application of diffuse-interface models for both single- 
component and binary fluids. These models have foundations in statistical mechanics, ki- 
netic theory, mechanical theories, and nonequilibrium thermodynamics. They provide an 
alternative approach to the classical hydrodynamic free-boundary problem and have been 
used successfully in many applications. They should continue to play an important role, in 
concert with other theoretical and experimental efforts, in the understanding of hydrody- 
namic phenomena associated with critical fluids and other flows involving complex interface 
morphologies and topological changes. 
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Figure 1: Moving contact line: This figure shows the constant density contours 
in the inner region. The frame of reference is fixed on the contact line region 
so that the bottom plate moves from left to right at an imposed velocity. In 
this case the imposed static contact angle at the wall is approximately 54° 
while the dynamic contact angle (away from the immediate vicinity of the 
contact line) is approximately 125°. The parameter values consistent with the 
notation given in Seppecher (1996) are iZ = 20, Um = 2, if = 10, y = -0.3, 
and Ca = 40 x 10~ 3 . (This figure was provided by P. Seppecher). 
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Figure 2: Moving contact line: This figure shows the streamlines in the inner 
region associated with the density contours shown in Fig. 1. The frame of 
reference is fixed on the contact line region so that the bottom plate moves from 
left to right at an imposed velocity. (This figure was provided by P. Seppecher). 
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Figure 3: Spinodal decomposition: This time sequence (top to bottom) shows 
spinod&l decomposition of a binary fluid in a rectangular cell across which is 
imposed a vertical temperature gradient (the cell is hottest at the bottom). 
The composition is indicated in greyscale. These computations show the sys- 
tem evolving from an initially random composition distribution in the top 
frame, to one in the bottom frame where the domain structure size varies with 
the temperature. (This figure was provided by D. Jasnow and J. Vinals). 
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Figure 4: Flow in a driven cavity: This time sequence (across the top, left 
to right and across the bottom, left to right) shows the order parameter (in 
greyscale) representing the binary fluid composition. The flow is set up by 
an imposed velocity (left to right) on the bottom surface of the cavity. The 
densities and viscosities of both phases are equal. (This figure was provided 
by R. Chella and J. Vinals). 
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